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1. 


1 .  INTRODUCTION 

This  paper  is  concerned  with  a  theoretical  and  numerical  investigation  of  the 
propagation  of  stress  waves  induced  by  surface  impact  on  a  laminated  plate  of  finite 
depth  and  infinite  lateral  extent.  Attention  is  focussed  on  the  problem  of  an  impulsive 
line  load  acting  on  the  upper  surface  of  the  plate  and  generating  plane  wave  disturbances 
travelling  in  the  laminate  along  the  direction  normal  to  the  line  load.  The  plate  is 
constructed  of  an  arrangement  of  layers  (or  plies)  of  fibre-reinforced  material  in  which 
the  reinforcement  of  each  layer  is  a  family  of  parallel  fibres  lying  in  the  plane  of  the 
layer.  The  plies  are  assembled  in  a  periodically  repeating  configuration  of  N  unit  cells, 
the  configuration  being  defined  in  terms  of  the  angles  between  the  fibre  direction  in  each 
ply  and  some  specified  reference  direction.  For  the  present  we  have  concerned  ourselves 
with  a  simple  0/90  configuration.  Here  the  unit  cell  consists  of  two  inner  layers,  each 
of  thickness  h,  with  the  fibres  running  along  the  reference  direction,  and  bounded  above 
and  below  by  a  layer  of  thickness  h  of  the  same  material  with  fibre  direction  orthogonal 
to  the  reference  direction  (see  Figure  17).  This  choice  of  unit  cell  is  made  purely  for 
simplicity  and  the  techniques  that  we  use  are  applicable  to  any  configuration  of  the  unit 
cell. 

We  model  the  fibre  reinforced  material  as  a  homogeneous  continuum  of  transversely 
isotropic  elastic  material  with  the  axis  of  transverse  isotropy  along  the  fibre  direction. 

This  means  that  we  look  at  waves  whose  wavelengths  are  an  order  of  magnitude  greater 
than  the  fibre  diameter  and  inter-fibre  spacing  so  that  on  the  scale  of  the  wavelength 
the  continuum  theory  might  be  expected  to  be  valid.  Typical  of  the  materials  we  have 
in  mind  is  the  ICI  product  PEEK,  formed  of  carbon  fibres  embedded  in  a  thermoplastic 
resin,  for  which  typical  dimensions  are  h  s  125  jun  with  the  fibre  diameter  and 
inter-fibre  spacing  of  the  order  of  6  nm.  Thus  we  are  thinking  in  terms  of  wavelengths 
of  the  order  of  1/2  to  1/3  the  ply  thickness  or  greater  for  which  the  non-dimensional 
wave  number  kh  =  2trh/A  (where  A  is  the  wavelength)  varies  from  zero  at  infinite 
wavelength  to  a  value  of  approximately  18  at  A  =  h/3.  For  smaller  wavelengths,  of  the 
order  of  h/10  or  less,  the  continuum  model  will  break  down  due  to  diffraction  and 
scattering  by  the  individual  fibres. 
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There  is  a  considerable  simplification  in  the  mathematics  to  be  gained  by  treating 
the  composite  as  inextensible  in  the  direction  of  transverse  isotropy.  This  is  an 
idealization  of  the  fact  that  the  extensional  modulus  of  the  continuum  along  the  fibre 
direction  can  be  of  the  order  of  100  times  that  in  the  cross  fibre  direction. 
Mathematically,  the  effect  of  the  idealization  is  to  reduce  the  order  of  the  differential 
equations  and  this  leads  to  solutions  involving  fewer  parameters.  A  consequence  of  this 
reduction  in  the  order  of  the  equations,  however,  is  that  it  is  no  longer  possible  to 
satisfy  all  the  interface  continuity  conditions  between  the  plies.  This  leads  to  a  singular 
perturbation  problem,  in  which  it  is  necessary  to  allow  the  tangential  component  of 
traction  along  the  fibre  direction  to  be  discontinuous  across  the  interface,  with  a 
consequent  singularity  in  the  stress  component  along  the  fibres,  associated  with  a  finite 
load  carried  by  the  surface  layer  of  fibres.  This  singular  perturbation  problem  has  been 
examined  in  detail  for  static  problems  by  Everstine  and  Pipkin  [1  ]  and  for  dynamic 
problems  by  Green  [2],  Green  and  Milosavljevif  [3]  and  Bayiis  and  Green  [4], [5].  In 
particular,  Bayiis  and  Green  [4], [5],  dealing  with  a  four  ply  laminate,  derive  an 
expression  relating  the  singularity  along  the  fibre  direction  to  the  discontinuity  in  shear 
stress.  They  present  detailed  comparisons  of  the  stresses  in  the  inextensible  laminate 
with  those  in  a  strongly  anisotropic  but  not  inextensible  material.  These  comparisons 
show  that  the  inextensible  theory  provides  an  acceptable  approximation  to  the  stress 
variation  through  the  laminate  except  in  the  very  long  wavelength  region  and  the 
dispersion  curves  for  the  two  material  models  are  also  shown  to  be  in  good  agreement, 
again  with  the  exception  of  the  very  long  wavelength  limit.  As  with  singular 
perturbation  problems  generally,  the  shear  stress  discontinuities  are  to  be  interpreted  in 
terms  of  very  narrow  bands  (boundary  layers)  adjacent  to  the  interfaces,  through  which 
there  exist  high  stress  gradients,  giving  large  changes  in  stress  across  the  bands.  The 
associated  singular  stresses  along  the  boundary  fibres  are  to  be  interpreted  as  high  stress 
levels  in  the  boundary  layers,  which  contribute  finite  loads  in  the  fibre  directions  when 
integrated  through  the  boundary  layers.  It  is  with  these  interpretations  in  mind,  that  we 
adopt  the  idealization  of  inextensibility  to  give  a  mathematically  simple  model  of  our 
fibre  reinforced  material.  Having  employed  this  simplified  model  to  develop  and  apply 
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our  analytical  and  numerical  techniques,  it  is  proposed  in  future  work  to  drop  the 
constraint  of  inextensibility  and  to  employ  the  full  system  of  governing  equations 
appropriate  to  the  continuum  model  of  the  composite  as  a  transversely  isotropic  elastic 
material. 

In  the  work  reported  here,  we  make  no  assumptions  about  the  variation  of 
displacements  and  stresses  through  the  laminate,  such  as  is  done  in  engineering  theories 
of  plates  and  shells.  Our  method  is  to  solve  exactly  the  system  of  governing  equations 
appropriate  to  each  layer,  matching  the  solutions  across  the  interfaces  and  satisfying  the 
appropriate  boundary  conditions  at  the  upper  and  lower  surfaces  of  the  laminate.  This 
analytical  solution  is  carried  out  in  Section  ?•  where  we  work  in  terms  of  systems  of 
coupled  ordinary  differential  equations,  obtained  by  taking  Laplace  transforms  in  time  and 
Fourier  transforms  in  the  in-plane  spatial  coordinates,  of  the  governing  equations  of  our 
model.  These  ordinary  differential  equations,  involving  the  coordinate  normal  to  the  plies 
as  independent  variable,  are  solved  using  the  propagator  method  of  Gilbert  and  Backus 
[6],  as  applied  to  the  inextensible  laminate  by  Green  &  Baylis  (7].  This  yields  the 
exact  solution  for  the  variation  of  the  transforms  with  depth  throughout  the  laminate  and 
.he  approximations  arise  only  in  the  numerical  methods  for  inverting  the  transforms. 

The  problem  of  inverting  the  transforms  is  considered  in  Section  3.  There  we  show 
that  it  is  necessary  to  solve  the  dispersion  equation  (relating  frequency  to  wavelength)  for 
plane  waves  travelling  in  the  plate  under  traction  free  boundary  conditions  at  the  upper 
and  lower  surfaces.  Inversion  also  requires  that  we  evaluate  the  residue  of  the  inversion 
integral  along  each  branch  of  the  dispersion  equation  and  we  must  perform  an  integration 
along  each  branch  and  a  summation  over  all  branches  to  obtain  the  solution.  In  this 
Section  we  also  briefly  outline  the  stationary  phase  approximation,  which  gives  the 
solution  at  large  times.  This  brings  out  the  role  of  the  group  velocity  in  the  asymptotic 
solution  and  draws  attention  to  the  wave  fronts,  which  are  associated  with  maxima  and 
minima  of  the  group  velocity  and  provide  the  dominant  contributions  at  large  times. 

The  numerical  solution  of  the  dispersion  equation  and  the  subsequent  numerical 
integration  over  the  branches  of  the  dispersion  curves  is  outlined  in  Section  4.  It  is 
here  that  the  approximations  come  into  play,  since  the  solutions  are  necessarily  limited  to 
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a  finite  number  of  branches  and  the  integrations  over  each  branch  must  be  limited  to  a 
finite  range  of  values  of  the  wave  number,  It.  We  have  chosen  as  an  upper  value 
k  =  20,  corresponding  to  wavelength  of  the  order  of  1/3  the  ply  thickness  but  6  times 
the  fibre  diameter  and  inter-fibre  spacing,  which  we  estimate  as  being  the  limit  at  which 
the  continuum  model  would  be  valid.  We  show  that  restricting  the  integral  to  a  finite 
range  of  values  of  k  gives  rise  to  the  phenomenon  of  "windowing"  and  we  apply  the 
technique  due  to  Hamming  [8]  in  an  attempt  to  reduce  this  effect. 

Section  5,  which  completes  the  paper,  contains  detailed  results  in  graphical  form  and 
some  discussion  of  the  contribution  of  the  higher  harmonics  as  compared  with  the  first 
two  modes.  The  results  consist  of  graphs  of  the  dispersion  equations,  which  illustrate  the 
plateau  and  step  regions  arising  in  the  higher  harmonics  and  of  the  corresponding  group 
velocity  curves,  which  display  the  contributions  of  the  high  harmonics  to  the  wave  front 
disturbance  at  large  times.  Graphs  are  also  presented  which  show  the  variation  of 
displacements  and  stresses,  at  the  upper  and  lower  surfaces  and  on  the  middle  plane,  as 
functions  of  position  at  various  times. 


2.  GOVERNING  EQUATIONS  AND  PROPAGATOR  SOLUTIONS 

We  consider  a  laminated  plate  composed  of  layers  of  the  idealized  fibre  reinforced 
material  which  is  inextensible  in  the  fibre  direction,  the  fibres  being  in  the  plane  of  the 
layers  and  with  the  fibre  directions  in  adjacent  layers  being  orthogonal  to  each  other. 

We  choose  a  Cartesian  Coordinate  system  of  axes  with  origin  in  the  middle  surface  of 
one  of  the  layers  and  with  Ox3  parallel  to  the  fibre  direction  in  that  layer.  Letting 
Ox2  be  parallel  to  the  fibre  direction  in  the  adjacent  layers  then  the  normal  to  the 
laminate  coincides  with  the  x ,  -axis.  We  shall  be  concerned  with  plane  waves  generated 
by  a  line  load  acting  on  the  upper  surface  of  the  laminate  along  a  line  making  an  angle 
(x/2— y)  with  the  x3-axis.  The  resulting  disturbance  will  generate  a  displacement  vector 
u(x ,  ,x2,x3,t)  of  the  form 

u(x,,x2,x3,t)  =  u(x,,x,t)  ,  (1) 

where 

x  =  x3  cos  y  +  x2  sin  7  ,  (2) 

which  corresponds  to  a  plane  wave  whose  normal  I  ts  in  the  x2x3  plane  at  an  angle  — y 
to  the  x3-axis  (see  Figure  18).  The  components  tjj(x,  ,x2,x3,t)  of  stress  in  each  layer 
of  the  laminate  will  likewise  reduce  to  functions  tjj(x,,x,t)  and  these  are  given  in  terms 
of  the  strain  components  via  the  constitutive  equations  appropriate  to  each  layer  (see  for 
example  Baylis  &  Green  [4]).  Before  proceeding  further  it  is  convenient  to  take  the 
Laplace  Transform  in  time  and  Fourier  Transform  with  respect  to  the  space  variable  x  of 
both  the  displacement  components  Uj(x,,x,t)  and  the  stress  components  tjj(x,,x,t).  Thus 
we  define  dispacement  transforms  U.V.W  by 

U(x,,k,s)  -  |  ^  |  u, (x, ,x, t )e~ste-ikx  dt  dx  , 

V(x,,k,s)  -  [*”  f  u2(x,  .x.tje'^e"'^  dt  dx  ,  (3) 

J  -00  J  0 

W(x.,k,s)  -  f  f  u, (x, , x, t je-*1 e~ ' dt  dx , 

J-oo  J0 

and  stress  transforms  Tjj(x,,k,s)  by 


and  the  stress  transforms  are  given  by 


(7) 


T22  -  p(c J-2c2)  ^  T,  +  iksp4c2[l  -  £f]v,  , 

1  1 

T23  -  pc^ikc  V,  ,  T31  -  pcMkc  U, 

T33  -  -P  ^  <T.  +  2iks  V,)  ■ 

1 

In  equations  (5)  and  (7),  p  is  the  density  of  the  continuum  and  c*,c£,c£  are  squared 
wave  speeds  derived  from  the  elastic  constants  of  the  continuum  (Green  [2]). 

The  solution  of  the  differential  equations  (4)  may  be  written  in  terms  of  the 
propagator  matrix  P(x,-x,)  in  the  form 

Y(x,)  =  P(x,-xt)Y(x,)  ,  (8) 

where  xt  is  some  fixed  value  of  x,  within  the  layer  and  x,  is  any  other  value  in  the 
same  layer,  Y(x,)  is  the  vector  (T^x,),  S,(x,),  U,(x,),  V^x,))^  where  T  denotes  the 
transpose,  and  the  components  of  the  4x4  matrix  P  are  given  in  Table  1,  page  33. 

Using  subscript  2  to  denote  the  transformed  quantities  within  a  layer  of  material 
with  the  fibre  direction  parallel  to  the  x2-axis,  the  constraint  condition  gives 

V2  =  0  ,  (9) 

and  the  equations  equivalent  to  (5)  involve  U2,  W2,  T2  and  R2  where  R2  =  T,3/pc'. 
Letting  Z(x,)  denote  the  vector  (T  2,R2,U  2,W  2)T  the  solution  of  these  equations  is  given 
by 

?(*,)  =  Q(x,  -*  ,)?(*,)  (1C) 

where  i,  is  some  fixed  value  of  x,  within  the  layer  and  x,  is  any  other  value  in  the 
same  layer.  The  components  of  the  4x4  propagator  matrix  Q  may  be  derived  from 
those  of  P  by  the  substitutions  given  in  Table  2,  page  33. 

Equations  (8)  and  (10)  yield  the  solution  of  the  governing  equations  in  individual 
layers  of  material  with  fibres  parallel  to  the  x3-axis  and  the  x2~axis  respectively.  In 
order  to  obtain  the  solution  in  a  multi-layered  laminate  it  is  necessary  to  satisfy  the 
appropriate  conditions  at  the  upper  and  lower  surfaces  as  well  as  continuity  conditions  at 
the  interfaces.  In  general,  the  conditions  at  the  interface  between  two  dissimilar  elastic 
materials  which  are  bonded  together  requires  continuity  of  all  three  displacement 
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components  and  the  three  components  of  traction  across  the  interface.  In  terms  of  the 
transformed  quantities  these  continuity  conditions  are 

U,  =  L'j.  V,  =  V2.  W,  =  W;. 

(11) 

T  ,  =  T  2  S ,  -  S  j.  a,  =  R  , 

For  the  idealized  inertensible  material,  however,  it  has  been  shown  by  Baylis  and  Green 
[4 1  that  there  exisita  the  possibility  of  a  discontinuity  in  the  tangential  component  of 
stress  parallel  to  the  fibre  direction  across  any  surface  in  the  material  and  in  particular 
this  allows  a  discontinuity  n  x  n  material  at  jny  interface  or  free  boundary  in  these 
materials  Thus  it  in  menace  between  materials  l  and  conditions  (11)  reduce  to 

i.  =  i.  r  .  :  v  i  vy  -  ■)  (U) 


the  last  two  conditions  nrinj  ■  .  nsrquencr  it  the  ineytensibilitv  constraints  applied  to 

the  dispiacemem  ontinuitv  ■  <n  .,ns 

i  .insider  now  i  .avei  t  tncinru  :ti  >t  material  /  sandwiched  between  two  layers 
it  material  .  citing  >  ■  mj  s  •  lenote  the  -aiue  'I  the  vector  Y  at  the  bottom  and 
top  of  this  avei  'espeitiveiv  nr  na.e  'tom  equation  Hi  that  * 

>'u  =  P'^hiY1  i  P'h.Wid'  (13) 

where  we  have  made  use  it  the  ptopertv  that  fot  anv  propagator  matrix 
P(4*-b|  =  PialPthiisee  Gilben  A  Backus  |b|)  If  we  now  impose  the  condition  that 
vy  s  0,  V*-  =  0.  arising  from  the  interface  conditions  (12).  n  is  possible  to  eliminate 
the  (discontinuous)  stresses  Su  S*-  from  equations  (13)  to  obtain 


Xu  -  XL 

det  R 


where 


r, 

r  ,  ,  r,;  ' 

rn  ri  i 

R  - 

|  ,  R  - 

u, 

r  i,  r22 

r2,  r, , 

(15) 


and  the  components  r,j  of  R  are 

given  in  terms 

of  the 

expressions 

ri  i  -  p,i  -  . 

ri  2  "  Pi  3  * 

pi  2S43 

Paa 

P42 

«■*,  -  P  3  1  -  • 

r  4  4 

r  2  2  -  P33  - 

Z22EA.3 

P42 

(16) 


9. 


In  the  same  way,  if  we  consider  a  layer  of  thickness  2h  of  material  2  sandwiched 
between  layers  of  material  /  we  may  express  the  value  of  X,  at  the  top  of  the  layer  in 
terms  of  the  value  at  the  bottom  in  the  form 


Applying  equation  (21)  to  a  sequence  of  (n-1)  unit  cells  then  gives 


10. 


x(rt-i)  _  0(n-l)  x(i) 


(\r‘-l>"'1)DX(i)  -  x1\,(x"-2-x"~2)x(^) 

(X, -X2) 

where  we  have  used  the  Cayley-Hamilton  theorem  to  express  Dn-1  in  terms  of  D,  I 
and  the  eigenvalues  X,,  X2  of  D.  Equation  (23)  refers  to  a  set  of  (n-1)  unit  cells 
assumed  to  be  embedded  in  a  repeating  pattern  so  that  the  constraint  conditions  arising 
from  inextensibility  is  operative  on  every  layer  of  thickness  2h  of  both  material  1  and 
material  2.  In  order  to  consider  wave  motion  in  a  laminate  of  finite  depth  2nh,  it  is 
necessary  to  encase  the  (n-1)  unit  cells  in  two  half  cells,  each  consisting  of  an  inner 
layer  of  material  2  and  an  outer  layer  of  material  1.  Then  each  of  the  layers  of 
material  2  is  still  constrained  by  the  inextensibility  condition  imposed  by  the  outer  layer 
of  material  /  and  the  expressions  previously  derived  for  the  transmission  through 
material  2  may  be  employed  at  the  top  and  the  bottom  of  the  (n-1)  cells  to  give 

x(n-i)  _  (f  pn-1  -  _L_jx(i>  .  ( 

Each  of  the  outer  layers  (of  thickness  h)  of  material  I  at  the  top  and  bottom  of  the 
complete  plate  must  be  treated  separately.  Each  is  subject  to  the  constraint  V,  =  0  at 
the  interface  with  material  2  and  it  will  be  assumed  that  the  tangential  component  of 
traction  S,  vanishes  at  the  outer  surface.  Applying  equation  (8)  to  the  layer  on  the 
upper  surface  of  the  plate  gives 

Y(h)  =  P(h)Y(0),  C 

and  making  use  of  the  conditions  S,(h)  =  0  and  V,(0)  =  0  allows  the  values  of  T,(h), 
U,(h)  to  be  given  in  terms  of  T,(0),  U,(0)  by  the  expression 

X,(h)  =  SX,(0)  ,  C 

where  the  elements  of  the  2  x  2  matrix  S  are  given  in  terms  of  the  elements  pjj  of 
P(h)  by 


S12  ~  Pi 


P] jPj2 
Pt, 


EiiEii 

Pl2 


l  2  -  P3  3 


.  P.3  zP  2  3 


(27) 


Applying  the  same  procedure  to  the  layer  on  the  bottom  surface  of  the  plate  leads  to 
the  expression 


X,(h)  -  S  X, (0)  -  ggf-g-  X, (0) 


and  the  last  expression  in  equation  (28)  uses  the  result  that  det  S  =  1 .  Combining  (24) 
with  (26)  and  (28)  gives  for  the  overall  laminate 

X(n)  -  S  F  Dn_  1  -r £-|  -  X(°) 

~  ~  ~  ~  det  F  S  “ 

-  M  X(0)  C 

where  x(®)  and  x(n)  denote  the  value  of  X  at  the  bottom  surface  and  the  top  surface 
respectively  and  M  is  the  overall  transmission  matrix  for  the  laminate.  If  the  upper 
surface  of  the  laminate  is  subjected  to  a  time  dependent  normal  line  load  of  the  form 
t,,(x2,x3,t)  =  P(t)  5(x3  cos  7  +  x,  sin  y)  =  P(t)  6(x) 
where  <5(x)  is  the  Dirac  delta  function,  and  the  lower  surface  of  the  laminate  is  traction 
free,  then 


X<n)  - 


X(0)  - 


where  F(s)  is  the  Laplace  transform  of  P(t).  Substituting  from  (31)  into  (30)  then  gives 
the  displacement  transform  at  the  bottom  surface  in  terms  of  the  stress  transform 


u(0)  .  ,  ( 

1  m12(k,s) 

and  the  displacement  transform  U(n)  at  the  bottom  surface  is  given  by 

l 

U(n)  _  p(§)  .  ( 

i  m, j (k, s) 

Equations  (32)  and  (33)  relate  to  the  normal  line  load  acting  on  the  upper  surface 
of  the  plate.  We  may  examine  the  case  of  a  tangential  line  load  in  a  similar  fashion. 
Because  of  the  inextensibility  constraint  in  the  x3-direction  at  the  upper  surface  the 
component  t ,  3  of  the  tangential  stress  in  the  x  3-direction  induces  a  singularity  in  the 


12. 


stress  component  t33  at  the  surface  but  transmits  no  disturbance  into  the  material.  It  is 
only  the  tangential  component  t12  in  the  x2 -direction  that  produces  a  wave  motion  in 
the  underlying  laminate.  To  obtain  the  solution  to  this  problem  we  return  to 
equation  (24).  This  relates  the  vector  x(i)  at  the  interface  between  the  bottom  layer  of 
material  /  and  the  rest  of  the  laminate  to  the  vector  x(n-i)  at  the  interface  between 
the  laminate  and  the  top  layer  of  material  I.  To  complete  the  solution  we  now  use 
equation  (25)  as  before  but  the  conditions  imposed  on  each  of  the  outer  layers  of 
material  1  are  that  V,  =  0  at  the  interface  of  the  layer  with  material  2  and  that 
T ,  =  0  at  the  outer  boundary.  This  allows  us  to  express  the  vector  Xu  =  (S ,  U ,  )T  at 
the  upper  surface  of  the  top  layer  in  terms  of  the  vector  X^-  =  (T ,  U ,  )T  between  the 
top  layer  and  the  rest  of  the  plate  in  the  form 

X  =  G  XL  ,  (34) 


where  the  elements  of  the  matrix  G  are  given  by 


8,i  -  Pa,  '  EliEzZ 
Pi  2 


P  2  3 


s2 


-  PllP-32 


822  -  P3 


-  P.IJ.P.22 
P>2 

_  Pi iPi 7 


(35) 


*31  n  on  fij  n 

Pi  2  P, 2 

Combining  equation  (34)  and  its  equivalent  at  the  bottom  layer  of  the  laminate,  with 
equation  (24)  then  yields 


x<")  - 

lr  v  nn-1  _ 

F 

G 

— 1  x(0) 

(36) 

where 

~  air 

F 

cJ  -  • 

c  - 

822  812' 

(37) 

822  8 i , 

and 

det  G  =  g, 

182  2  -  82181 

2 

= 

-  P 2  27P  1  2- 

(38) 

In  terms  of  the  overall  matrix  M  defined  by  equation  (30)  we  may  write  equation  (36) 


in  the  form 


X(n)  -  [c  S-'  M  |e’  |]  X<°>  -  N  X<°>  .  (39) 

For  a  tangential  line  load  with  component  t ,  t  given  by 

t,2(x3,x3.t)  =  Q(t)  S(x3  cos  7  +  x2  sin  7)  =  Q(t)  6(x)  (40) 

acting  on  the  upper  surface  of  the  plate  and  with  the  lower  surface  traction  free, 
equation  (39)  gives 


13. 


u(0)  .  _Q(|> 
i  n,  (k,s) 

(4 

0(n)  _  n22(k,S)Q(s) 

'  n, 2(k,s) 

where  5(S)  is  the  Laplace  transform  of  Q(t)  and  uW,  are  the  transforms  of  the 
normal  displacements  on  the  lower  and  upper  surfaces  respectively.  Using  the 
expressions  for  the  components  of  G  and  S  given  in  equation  (35)  and  (27)  the 
components  n12  and  n22  of  the  matrix  N  may  be  expressed  in  terms  of  the  components 


of  M  and  the  propagator  P(h)  in  the  form 


_  P32P1 2 


Combining  equations  (41)  linearly  with  (32)  and  (33)  gives  the  solution  for  an  arbitrary 
line  load  on  the  upper  surface  of  the  plate. 

Equation  (33)  gives  the  transform  of  the  normal  component  of  displacement  at  the 
upper  surface  due  to  the  specified  normal  line  load.  From  the  propagator  solution  (8) 
and  equations  (7)  it  is  possible  to  obtain  expressions  for  the  transform  of  the  tangential 
displacement  v(n)  at  the  upper  surface  and  the  transformed  components  of  stress.  Thes' 
are  given  by 


P«.p(§> 


(p„4  *  0) 


AV  -  P^{(l  -  %]  PCs)  +  iks  [l  -  £*-]  V<n)}  , 


-pc, IkcV, 


-r(n)  J,.  „(n) 

T;3'  -  pc3ikcU, 


+  2 iksV 


in)]  -  u'n)5(x,-2nh)} 


In  equation  (43)  the  component  T(n)  relates  to  the  limiting  value  on  approaching  the 

1 1 

surface  from  within  the  material  and  this  stress  component  jumps  discontinuously  to  zero 
on  crossing  the  surface.  The  Dirac  delta  function  term  in  the  expression  for  T'n)  gives 


the  singularity  in  the  reaction  stress  along  the  fibres  in  the  surface,  which  is  required  in 
order  to  balance  the  shear  stress  discontinuity.  (In  the  argument  of  the  delta  function, 


the  origin  of  coordinates  has  been  taken  to  be  in  the  middle  surface  of  the  2n  ply 
laminate.)  Results  appropriate  to  the  lower  surface  of  the  laminate  may  be  derived  from 


r(0 


r<*> 


PC2 


T<*3  -0-  T 


r(0 
1  3  3 


2  2 
P<-'3C  2 


T( 1 ) _T^  2 ) 


)  5(X*) , 

(48) 


2  2 
PC.3C2 


'.-F 


(t)_T(2) 


iks 


]«(**>, 


r<2> 


0, 


r(2) 


-  pc 


![■-%] 


where  the  argument  of  the  S  function  is  x*  =  x-(2n-l)h,  corresponding  to  the  origin  at 
the  middle  surface  of  the  laminate.  Results  similar  to  those  detailed  in  equations  (45)  - 
(48)  may  be  obtained  at  any  other  interface  of  the  laminate  whilst  at  any  interior  point 
in  any  layer  it  is  possible  to  also  determine  the  non-zero  in-plane  displacement  (V1  or 
W2).  In  the  next  section  we  consider  the  problem  of  inverting  these  transforms  to  give 
the  required  solutions. 


3.  TRANSFORM  INVERSION 


I 

I 


The  techniques  developed  in  Section  2  yield  the  transforms  of  the  displacements  and 
stresses  throughout  the  laminate  as  known  functions  of  the  transform  parameters  k  and  s 
at  any  value  of  xr  In  order  to  determine  the  displacements  and  stresses  as  functions  of 
the  coordinate  x  =  x3  cos  7  +  x2  sin  7,  normal  to  the  line  load  and  of  the  time  t,  it 
is  necessary  to  invert  these  transforms.  Typical  of  the  quantities  to  be  considered  is  the 
transform  of  the  normal  displacement  at  the  upper  surface  due  to  a  normal  line  load, 
which  is  given  by  equation  (33)  in  the  form 


0(n)  _  m2 •  f >  p(g) 
1  m,2(k,s) 


(33) 


Letting  u,(2nh,x,t)  denote  the  normal  displacement  on  the  upper  surface  x,  =  2nh, 
due  to  a  line  load  which  consists  of  a  unit  delta  function  in  time  P(t)  =  6(t),  for  which 
P(s)  =  1,  the  displacement  u^(x,t)  corresponding  to  the  transform  (33)  associated  with 
any  P(t)  is  then  given  as  the  convolution  of  P(t)  with  u,(2nh,x,t)  in  the  form 
P  f 1 

u, (x, t )  -  J  u, (2nh,x,r)P(t-r)  dr.  (49) 

Accordingly  we  restrict  attention  to  inverting  the  transform  (33)  with  P(s)  =  1 ,  for  which 
the  formal  solution  is 


,,  .  1  f”  p'+io°  m, . 

u.(2nh,x,t)  -  ■  1 .  , .  — 

i  ■*  -00  y—  ioa  nt* 


(KA1  est  e 


ikx 


ds  dk  . 


(50) 


'7-  too  1  2  (k  .  S) 

The  integral  with  respect  to  s  may  be  evaluated  in  terms  of  the  residues  of  the 
integrand  at  the  zeros  of  the  function  m ,  2(k,s)  in  the  left  half  plane.  The  equation 

m,2(k,iu)  -  0  (51) 

is  the  dispersion  equation  for  plane  wave  propagation  in  the  laminate,  corresponding  to 
waves  travelling  in  the  direction  of  the  normal  to  the  line  load  under  traction  free 
conditions  at  the  two  surfaces  of  the  plate.  This  equation  has  an  infinite  number  of 
pairs  of  roots,  ojj  =  twj(k)  (j=l,2,...),  each  pair  corresponding  to  forward  and  backward 
travelling  waves  associated  with  one  branch  of  the  dispersion  curve.  In  terms  of  these 
solutions,  equation  (49)  becomes 


u,  (2nh 


,x,t)  -  i-  f  dk  l  eKkx-istl 

ldm12/ds  I.  _  stuj 


(52) 


(k). 
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Both  m12(k,S)  and  m  2  2(k,s)  are  even  /unctions  of  $  and  equation  (52)  may  be  written  as 


where 


u.(2nh,x,t)  -  -  V  f  R;(k)  sin  u>;(k)t 
T  j-1  J  J 


,ikx 


dk 


?-(k)  -  L 

rK;  Ldm, 2/ds  L 


s  -  +Iu)i(k) 

1 


(53) 


(54) 


rm,,(k,s) 
Ldm, j/ds 


s  -  -iuij(k)  . 


It  may  also  be  shown  that  Rj(k)  is  an  even  function  of  k  and  equation  (53)  may  be 
further  simplified  to  give 


u, (2nh, x, t ) 


sin  oij(k)t  cos  kx  dk 


Equation  (43)  contains  an  expression  for  vOj)  which  is  valid  provided  p44  *  0  and 
which  therefore  has  the  same  singularities  as  U(n)  for  the  case  P(s)  =  1 .  It  may  be 
shown  that  p43/p44  is  an  odd  function  of  k  and  we  may  therefore  express  the  tangential 
displacement  v,(2nh,x,t)  associated  with  a  unit  delta  function  normal  line  load  as 


v,(2nh,x,t)  -  -  -  V  f  h  i(k)R  :(k)sin  cjj(k)t  sin  kx  dk  , 

T  Jq  J  J  J 

where  hj(k)  is  the  value  of  the  ratio  p43/p44  evaluated  at  s  =  iuj(k). 

The  formulae  for  the  transforms  of  the  stress  components  in  equation  (43)  involve  terms 
of  the  form  ikCK")  and  ikV^)  and  these  when  inverted  give  du,(2nh,x,t)/dx  and 
dv ,  (2nh,x,t)/dx  respectively. 

Equations  (45)  give  the  transformed  normal  displacement  and  normal  component  of 
traction  at  the  interface  between  the  top  surface  layer  of  material  I  and  the  rest  of  the 
laminate.  These  have  the  same  singularities  as  the  expression  (33)  and  their  inverses  yield 


r  i  2i 

u,  ^(2n-l  )h,  x,  t  J  -  — -  J  J  s,,(k)Rj(k)  sin  oij(k)t  cos  kx  dk, 

j-1  0 


(55) 


t, , [(2n-l )h , x, t] 


r“ 

I  Jg  s,2(k)Rj(k)  sin  ou  j  (kit  sin  kx  dk, 


where  we  have  used  the  fact  that  s,  ,(k)  is  an  even  function  of  k  and  s,2(k)  an  odd 

Equations  (46),  (47)  and  (48)  then  yield  the  remaining  stress  components 


function  of  k. 


on  the  two  sides  of  the  interface  in  a  similar  fashion,  the  resulting  integrals  being  of  the 
forms  given  in  equations  (55)  but  with  different  factors  multiplying  Rj(k). 

All  the  results  derived  in  this  section  relate  to  the  problem  of  the  normal  line  load 
acting  on  the  upper  surface  of  the  laminate.  Results  for  the  effect  of  the  tangential  line 
load  may  be  obtained  in  a  completely  analagous  way  by  starting  with  the  solutions  given 
in  equations  (41)  rather  than  using  the  solutions  (33).  Referring  to  equations  (42),  it 
may  be  seen  that  the  zeros  of  the  element  n12,  which  determine  the  singularities  of  the 
transforms,  occur  at  the  zeros  of  m)2  and  possibly  at  the  zeros  of  p22.  A  detailed 
examination  of  the  product  p22m)2  shows,  however,  that  this  does  not  vanish  at  the 
zeros  of  p  2  2  and  therefore  the  singularities  of  these  transforms  occur  on  the  same 
dispersion  curves  as  for  the  normal  line  load  problem.  Thus  all  the  inversion  integrals 
associated  with  the  tangential  line  load  problem  involve  the  residues  Rj(k)  multiplied  by 
some  appropriate  factor.  Hence  for  both  the  normal  and  tangential  line  loads  the 
problem  of  inverting  the  transforms  reduces  to  evaluating  an  infinite  sum  of  infinite 
integrals  of  the  form 

I  C  fj(k)Rj(k)  si"‘\i<k>t  Co"3dk  (56) 

J-l  u 

where  fj(k)  is  some  factor  arising  from  a  function  of  k  and  s  which  is  evaluated  on  the 
branch  s  =  iuij(k)  of  the  dispersion  curve. 

The  expression  (56)  consists  of  a  sum  of  integrals,  one  along  each  branch  of  the 
dispersion  curve.  In  general  both  the  integration  and  summation  have  to  be  carried  out 
numerically  and  we  must  therefore  limit  the  range  of  integration  to  some  finite  interval 
(0,k)  and  restrict  the  summation  to  a  finite  number  of  branches  j=l,...,P  of  the 
dispersion  curve.  We  then  have  to  construct  a  computer  programme  to  solve  the 
dispersion  equation 

m,  2(k,i<i)j)  =  0  (57) 

numerically  in  order  to  obtain  values  of  uj(k)  along  each  of  the  P  branches  for  values  of 
k  taken  at  M  intervals  Ak  =  ic/M  from  k=0  to  k=k.  It  is  also  necessary  to  evaluate  the 
residues  Rj(k)  at  each  of  these  values  of  k  for  each  of  the  branches  j=l,2,...P. 

An  alternative  approximate  solution,  valid  for  large  times,  may  be  derived  by  the 


method  of  stationary  phase.  Writing  the  integrals  which  appear  in  equation  (52)  in  the 
form 


f”  Rj(k)  e*  (kx-ujt)  dk  _  T  R :  (k)  e  lt9J  dk  , 

1  -00  J  '  -oo  J 

where 

9j  -  (Uj(k)  -  k  x/t) , 

the  rapid  oscillation  produced  by  the  exponential  term  for  small  changes  in  k  (and 
therefore  9j)  at  large  values  of  t  produces  cancellation  of  the  integrands,  except  for 
those  values  of  k  for  which  dtfj/dk  is  zero.  The  major  contribution  to  the  integral  then 
comes  at  those  values  kg  of  k  for  which  is  stationary  and  which  are  given  by  solving 
the  equation 


Cj (kg)  -  x/t  ,  (58) 

where  Cj(k)  =  dajj/dk  is  the  group  velocity  on  the  j111  branch  of  the  dispersion  curves. 

For  a  specified  value  of  t,  equation  (58)  determines  one  or  more  values  of  kg  at  any 
given  value  of  x  and  the  contribution  from  the  integral  at  that  value  of  x  and  t  is  then 
given  approximately  by  * 

I* 


2 x  i 
d2oi : 


dk 


Rj (kg)  e' (kg_uj(kg)t  , 


g 

provided  d  Jo)j/dkg  *  0.  At  the  stationary  points  of  the  group  velocity, 

dCj/dk  =  dJwj/dk*  =  0,  this  result  breaks  down  and  the  approximation  must  be  carried 

to  higher  order,  giving  the  result 


d3u) ; 
dvT’ 
g 


Rj (kg) 


*  i  (kgX-oij  (kg)  t ) 


(591 


Thus  the  stationary  phase  approximation  at  ordinary  points  of  the  group  velocity  curve 
shows  that  the  long  time  solution  deacays  as  whereas  at  maxima  and  minima  of  the 
group  velocity  the  decay  goes  more  slowly  as  t-'  /3.  These  are  the  points  corresponding 
to  the  wave  fronts  and  it  is  these  wave  front  solutions,  given  by  equation  (59),  which 
dominate  at  large  times.  A  detailed  discussion  of  these  phenomena  for  a  single  plate  of 
isotropic  elastic  material  is  given  in  the  paper  by  Jones  [9], 


4.  NUMERICAL  METHODS 


To  carry  out  the  numerical  evaluation  we  choose  the  values  for  the  material 
parameters  that  were  previously  employed  by  Green  &  Baylis  [4]  and  which  are  derived 
from  measurements  carried  out  by  Markham  [10]  on  a  carbon  fibre/epoxy  resin 

composite.  For  the  inextensible  model,  these  become  c2/c2  =  4.297  and 

1  2 

c2/c2  =  2.301.  We  must  also  specify  the  number  N,  of  unit  cells  forming  the  laminate 

3  2 

and  the  results  reported  here  relate  to  the  simplest  case  of  N=1  although  the  computer 
programmes  have  been  written  to  cope  with  the  general  case  of  any  specified  number  N 
of  unit  cells.  In  order  to  obtain  the  dispersion  equation  (57)  and  the  expression  for  the 
residues  given  in  equation  (54),  we  have  made  use  of  the  algebraic  manipulation 
programme  REDUCE.  We  provide  as  basic  inputs,  the  elements  p;j  of  the  propagator 
matrix  P(h)  defined  in  Table  1  and  the  elements  qrs  of  the  matrix  Q(h),  obtained  by 
the  transformations  given  in  Table  2.  The  programme  then  calculates  the  elements  of 
the  matrices  R,  R,  F,  F,  using  equations  (16)  for  rjj  and  equivalent  equations  for  f;j, 
and  compounds  these  to  produce  the  matrices  C,  C  and  D.  It  may  be  seen  from  the 
definition  of  D  given  in  equation  (22)  that  det  D  =  1  so  that  the  eigenvalues  X ,  and  X  2 
of  D  satisfy  the  relation  X ,  X  2  =  1 .  It  is  also  a  consequence  of  the  definition  that 
d2J=  dt,  and  making  use  of  these  two  results  yields  explicit  expressions  for  X,  and  X2 
as 


subroutine  to  evaluate  m,2(k,$)  and  a  subroutine  for  R(k,§)  =  m  2  2/(dm ,  2/df)  from  which 
to  calculate  the  residues  Rj(k)  and  these  subroutines  are  incorporated  into  the  computer 
programme  which  solves  the  dispersion  equation. 

The  computer  programme  to  produce  the  solutions  u)j(k)  to  the  dispersion  equation 
ml2(k,Sj(k)  =  0,  where  sj (k)  =  icjj(k),  is  based  on  the  notion  of  fixing  k,  and  then, 
starting  from  u>  =  0,  marching  up  the  oi-axis  evaluating  m,  2  at  each  step  until  the 
required  number  of  zeros  of  m,2  have  been  determined.  A  zero  is  indicated  by  a 
change  of  sign  of  m12  at  two  consecutive  values  of  u>,  and  by  reducing  the  step  length, 
this  zero  can  be  determined  to  any  given  degree  of  accuracy. 

Clearly,  this  method  is  open  to  error  in  that  two  changes  of  sign  which  occur 
within  the  same  step  length  will  be  missed  and  a  change  of  sign  indicates  the  presence 
of  an  odd  number  of  zeros,  not  necessarily  just  one.  These  two  possibilities  did  in  fact 
give  rise  to  problems  since  adjacent  harmonics  do,  on  occasion,  run  very  closely 
together.  In  theory,  the  answer  is  simple  -  the  step  length  has  to  be  chosen  sufficiently 
small!  However,  since  a  substantial  number  of  harmonics  are  required,  covering  a  very 
large  set  of  values  of  k,  the  programme  would  then  become  excessively  expensive  in 
computer  time. 

A  refinement  introduced  to  save  running  time  is  to  estimate  a  value  of  uij(k)  and 
then  proceed  to  a  fine  search  for  a  zero  in  its  locality.  The  programme  is  initialized  by 
using  the  value  for  uj(0),  derived  by  solving  the  dispersion  equation  at  k  =  0 
analytically.  It  is  then  possible  to  use  the  history  of  the  j**1  harmonic  to  estimate  the 
location  of  the  j1*1  root  under  inspection.  For  each  new  value  kj+j  of  k  the  root 
corresponding  to  the  fundamental  mode,  j  =  1,  is  the  first  to  be  located  and  this  is 
estimated  by  u  =  u),(kj)  -  e,  where  r  =  56k,  and  6k  is  the  step  length  in  k.  This 
allows  for  the  possibility  of  the  dispersion  curve  having  a  negative  gradient.  The  step 
length  h,(kj+j)  in  o>  is  then  taken  as  the  minimum  of  0.1  e  and  0.1  {cu 2(kj)  -  u>,(kj)}. 
Thus,  if  the  previous  history  indicates  that  the  fundamental  mode  and  first  harmonic  are 
very  close  at  kj,  the  step  length  at  k;+1  is  chosen  as  a  tenth  of  the  gap  at  k;. 

Otherwise,  experience  indicates  that  an  increment  of  e/10  is  sufficiently  small. 

For  each  of  the  remaining  roots,  uij(kj+j)  for  j  >  1,  the  first  estimate  is  based  not 
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only  on  the  history  of  that  particular  harmonic,  but  also  on  the  roots  already  located  at 
kj+[  and  is  given  by 

u  =  max{uj(kj)  -  e.  oj_t(ki+1)  +  0.1  hj_j(ki+i)}. 

That  is,  if  the  root  of  uij(k;)  is  sufficiently  far  above  the  previous  root  at  k;+j, 
Uj_i(kj+i),  a  substantial  amount  of  computer  time  can  be  saved  by  stepping  the  first 
estimate  of  u)j(kj+j)  over  this  gap.  The  step  length  for  locating  this  root  is  determined 
by 

hj(k1+1)  -  min  {«/10  ,  uJtl<!yj<kl>  ,  }. 


This  allows  for  the  possibility  of  the  j1*1  harmonic  at  kj  being  very  close  to  either  the 
(j+l)th  harmonic  or  the  (j-1)1*1  harmonic. 

Once  the  step  length  hj  and  a  first  estimate  ui  for  oij(kj+j)  have  been  established 
the  procedure  for  locating  the  root  «j(kj+j)  is  as  follows:  given  kj+i  and  ui,  m,  2  can 
be  evaluated,  to  is  then  increased  by  the  step  length  hj  and  m ,  2  recalculated  with  this 
new  value  of  u.  We  continue  to  increase  oj  by  the  step  length  until  a  change  of  sign 
of  m,j  is  observed.  The  step  length  is  then  reduced  by  a  factor  of  10,  and  the  process 
is  repeated  using  the  last  value  of  ui  before  the  change  of  sign  occurred  as  a  new 
estimate  of  the  root.  Once  a  change  of  sign  has  been  re-established,  the  step  length  is 
reduced  by  a  further  factor  of  ten  and  the  complete  procedure  repeated,  with  termination 
occurring  when  oij(k)  has  been  determined  to  the  given  degree  of  accuracy. 

Having  obtained  the  solutions  of  the  dispersion  equation  and  the  associated  residues, 
we  are  now  in  a  position  to  perform  the  numerical  integrations  and  summation  involved 
in  inverting  the  transforms  through  expressions  of  the  form  (58).  To  do  this,  we  restrict 
the  range  of  summation  from  j  =  1  to  P  and  the  interval  of  integration  from  k  =  0  to 
k  and  rewrite  the  expression  (58)  in  the  approximate  form 


2 

j-i 


|  fj(k)Rj(k)  sin  u>j(k)t 


cos  kxl 
sin  kxJ 


dk 


sin  uij (k) t 


cos  kxl 
sin  kxJ 


dk 


J0  H(k'°  { 


cos  kx 
sin  kx 


}  dk 


(62) 


where 
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P 

H(k, t )  -  2  fjOORjOO  sin  (jj(k)t .  (63) 

j-1 

Of  the  two  approximations  involved  in  equation  (62),  the  effect  of  restricting  the 
summation  to  a  finite  number  of  roots  (P)  cannot  be  evaluated  without  some  estimate  of 
the  contributions  arising  from  the  residues  of  the  excluded  branches  and  in  general  such 
an  estimate  is  not  available.  The  effect  of  restricting  the  integral  to  a  finite  range  is, 
however,  immediately  calculable.  To  do  this  it  is  convenient  to  define  the  function 
G(k,t)  such  that 

C  H(k'°  (sin  3  dk-  *  C.c<k*«>  eikX  dk  '  (64) 

so  that  if  the  integrand  is  H(k,t)cos  kx,  then  G(k,t)  =  H(k,t)  for  k  >  0  and 
G(k,t)  =  H(-k,t)  for  k  <  0,  whereas  if  the  integrand  is  H(k,t)sin  kx  then 
G(k,t)  =  H(k,t)  for  k  >  0  and  G(k,t)  =  -H(-k,t)  for  k  <  0.  Letting  F(x,t)  denote  the 
function  whose  Fourier  transform  is  G(k,t),  we  have  that 

C(k, t )  -  T  F(y, t)  eiky  dy  (65) 

J  -CO 

and  using  this  gives 


J0H<k’t>  Cin  3  OK-  *  J”  G(k.t)  eikx  dx 


i  [  F(y,t)  f  eik(x-y)  dk  dy 
—C&  J  r. 


11  r(y.«)  -y. 


-  F(x,t)  .  (66) 

The  function  F(x,t)  defined  in  equation  (66)  is  the  convolution  of  the  true  signal  F(x,t) 
with  the  function  (sin  kx)/x  and  exhibits  the  phenomenon  of  "windowing".  This 
phenomenon  produces  a  spurious  oscillation  of  wavelength  2i/k  and  in  order  to  reduce 
the  effect  of  the  oscillation  we  have  made  use  of  the  Hamming  window  function  s(k) 
defined  by 
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s(k)  -  a  +  (1-a)  cos  fj^  1  ,  (67) 

lk  1 

where  a  is  some  parameter  satisfying  O  <  a  <  1 .  The  procedure  is  to  replace  the 
integrand  over  the  finite  range  (0,k)  by  the  product  of  the  integrand  with  s(k)  to  give  a 
new  approximation  F^(x,t)  given  by 

F„(x,t)  -  J*  s(k)H(k,t)  dk  . 

k 

-  i  f  s(k)C(k,t)dk  , 

-k 

-  a  F(x,  t)  +  (f(x-  I,t)  +  F(x  +  J,t)]  .  (68) 

The  convolution  terms  in  the  braces  are  each  half  a  wavelength  out  of  phase  with  the 
convolution  F(x.t)  and  this  serves  to  dampen  out  the  oscillation.  In  evaluating  our 
numerical  results  we  have  used  the  value  a  =  0.54  recommended  by  Hamming  [8].  Note 
that  both  F(x,t)  and  Fpj(x.t)  tend  to  the  true  signal  as  k  -t  ®.  There  is  one  further 
approximation  required  in  order  to  perform  the  integration  and  this  consists  of  replacing 
the  integral  by  a  finite  sum  of  terms.  This  has  been  done  using  the  trapezium  rule, 
with  interval  length  k/M  and  it  can  be  shown  that  this  procedure  is  equivalent  to 
replacing  the  convolution  integral  F(x,t)  defined  by  equation  (66)  by  the  convolution 
F(x,t)  defined  by 


F(x,t) 


co  sin  k(x-y)cos  ^  (x-y) 

F(y,t)  - * -  dy 

*  —CD 


sin  ^  (x-y) 

The  calculated  result  is  then  the  Hamming  function  of  this  convolution  and  is  given  by 
FH(x,t)  -  aF(x ,  t )  +  i!y*2.  [f(x  -  jr.t)  +  F(x  +  J,t}  . 


(69) 


(70) 


The  convolution  F(x,t)  tends  to  F(x,t)  in  the  limit  as  M  increases  to  infinity,  but  for 
finite  values  of  M  the  function  defined  by  equation  (69)  is  periodic  in  x  with  wavelength 
4xM/k  and  the  numerical  integration  procedure  is  therefore  limited  to  values  of  x 
satisfying  0  <  x  <  4rM/k. 
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In  performing  the  numerical  integration  we  have  chosen  as  unit  of  length  the  half 
thickness  h  of  each  ply  and  as  unit  of  time  the  quantity  h/c,.  For  a  given  value  of  7, 
the  dispersion  equation  has  been  solved  for  eighteen  modes  (P=18)  with  values  of  kh 
ranging  from  zero  to  20  (k=20/h)  in  steps  of  0.002,  corresponding  to  M  =  10,000.  This 
gives  rise  to  180,000  values  of  uij(k)  and  an  equal  number  of  values  of  the  residues  Rj(k) 
which  form  the  data  matrix  for  the  numerical  inversions.  The  procedure  for  this  is  to 
specify  a  value  of  t  and  to  form  the  sums  H(k,t)  appropriate  to  the  particular 
displacement  or  stress  being  evaluated.  Each  sum  H(k,t)  is  multiplied  by  the 
corresponding  value  of  the  Hamming  factor  s(k)  and  the  product  is  stored.  Using  the 
trapezium  rule,  the  integral  is  then  evaluted  for  a  range  of  values  of  x  from  x  =  0  to 
x  =  c,t  in  steps  of  Ax  =  2x/k. 


i 


5.  RESULTS 


The  results  we  present  in  this  section  are  in  the  form  of  graphs  and  fall  into  three 
sets.  The  first  set  consists  of  the  dispersion  curves,  relating  the  scaled  frequency  ofvc  2 
to  the  non-dimensional  wave  number  kh  for  a  range  of  values  of  the  angle  of 
propagation  y  satisfying  0°  <  y  <  90°.  Propagation  at  angles  outside  this  range  may  be 
obtained  from  these  by  symmetry.  The  dispersion  curves  together  with  their  associated 
residues  provide  the  fundamental  information  required  to  invert  the  transforms,  but  we 
have  seen  in  Section  3  that  the  long  time  solution  is  governed  by  the  group  velocity 
curves  and  in  particular  that  the  wave  fronts  are  related  to  the  maxima  and  minima  of 
the  group  velocity.  Our  second  set  of  results  displays  a  selection  of  group  velocity 
curves  at  different  values  of  the  angle  of  propagation  y.  Finally  we  present  a  set  of 
curves  showing  the  variation  of  displacements  and  stresses  at  the  outer  surfaces  of  the 
plate  as  functions  of  the  propagation  distance  x  at  various  times.  These  results  are 
derived  by  numerical  integration  along  the  dispersion  curves  and  are  presented  here  for 
propagation  at  the  angle  y  -  60°  only.  Calculation  of  each  of  these  latter  sets  of 
results  involves  the  use  Of  360,000  stored  values  for  each  angle  of  propagation  y  and 
results  for  other  angles  of  propagation  will  appear  elsewhere. 

The  laminate  with  which  we  are  dealing  consists  of  an  arrangement  of  4  plies  which 
is  symmetric  about  the  middle  surface.  It  may  be  shown  that,  in  consequence,  the 
dispersion  equation  factorizes  into  two  distinct  equations,  one  associated  with  the 
symmetric  (longitudinal)  motion  of  the  plate  and  the  other  with  the  antisymmetric 
(flexural)  motion.  The  dispersion  curves  for  the  fundamental  modes  of  these  two 
equations  have  been  examined  in  detail  by  Baylis  [11]  and  Baylis  and  Green  [4], [5]. 

They  show  that  the  limiting  velocity  of  short  waves  (large  kh)  propagating  at  angle  y  is 
either  the  velocity  of  Rayleigh  type  surface  waves  at  angle  y  in  the  outer  material  or  of 
shear  waves  at  angle  y  in  the  inner  material  according  as  to  whether  y  is  less  than  or 
greater  than  some  critical  value  yc,  which,  for  the  parameters  employed  here,  is  46.3°. 
The  limiting  short  wavelength  velocity  of  all  the  other  harmonics  is  either  the  velocity  of 
shear  waves  at  angle  y  in  the  outer  material  or  the  velocity  of  shear  waves  at  an  angle 
y  in  the  inner  material,  according  as  to  whether  y  is  less  than  or  greater  than  45°. 
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These  are  the  speeds  vls  =  (c^c2+c^s2)i  in  the  inner  material  (material  /)  and 

v2J  =  (cls2-*:^2)*  in  the  outer  material  (material  2)  at  which  p2  =  0  and  q2  =  0 

respectively.  We  can  also  associate  with  each  material  a  dilatational  wave  speed 

v,d  =  (c2s2+c|c2)i  and  v2d  =  (c2c2+c2s2)i  at  which  p,  =  0  and  q,  =  0,  respectively. 

For  values  of  y  <  45°  we  have  that  vld  <  v2d  but  v,s  >  v25  with  the  inequalities 
reversed  for  7  >  45°.  (Note  that  v)d  >  vM  and  v2d  >  v2s  for  all  7.)  On  any  branch 
of  a  dispersion  curve  at  points  where  the  phase  velocity  v  =  oj/k  is  greater  than  the 
supremum  (vld,v2d)  ail  of  p , ,  p2,  qt.  q2  are  pure  imaginary  and  the  solutions 
correspond  to  progressive  waves  in  all  regions.  As  the  phase  velocity  drops  below  the 
supremum  (vld,  v2d)  either  p,  or  q,  becomes  real  and  the  associated  dilatational  wave 
motion  is  evanescent.  Continued  reduction  in  v  makes  both  p,  and  q,  real, 
corresponding  to  evanescent  dilatational  disturbances  in  both  materials  and  as  v  decreases 
further  with  increasing  k  the  supremum  (vM,  v2J)  is  reached  making  either  p2  or  q2 
real  and  the  motion  subequently  involves  a  progressing  shear  wave  in  one  material  with 
all  the  other  waves  being  evanescent. 

Figure  1  shows  a  total  of  26  branches  (13  each  of  the  symmetric  and  antisymmetric 
modes)  of  the  dispersion  curves  for  propagation  at  the  angle  7  =  90°.  All  modes  except 
the  fundamental  mode  of  the  antisymmetric  motion  exhibit  a  cut-off  frequency  in  the 
long  wavelength  limit  kh  -*  0.  A  striking  feature  of  the  graph  is  the  existence  of  the 
two  clear  ghost  lines  brought  about  by  the  osculation  of  the  branches.  (These  are 
particularly  clear  when  the  graph  is  viewed  at  almost  grazing  incidence).  These  ghost 
lines  have  slopes  oo/k  =  c,  and  oi/k  =  c3,  which  correspond  to  v  =  vld  and 
v  =  v2d  =  v2J  respectively  at  7  =  90°.  As  the  dispersion  curves  approach  the  first 
ghost  line  from  the  left  they  exhibit  the  plateau  and  step  phenomenon  described  by 
Redwood  [12].  Along  the  plateaux  the  curves  are  almost  parallel  to  the  ghost  line,  with 
the  phase  velocity  and  group  velocity  virtually  constant  at  the  value  given  by  the  slope 
of  the  ghost  line.  On  the  steps  the  phase  velocity  exhibits  a  small  but  sudden  drop, 
with  an  associated  drop  in  the  group  velocity,  either  ending  up  on  the  next  plateau,  if 
still  to  the  left  of  the  ghost  line,  or  passing  through  and  moving  rapidly  towards  the 
second  ghost  line,  where  the  phenomenon  is  repeated.  On  crossing  the  second  ghost  line 


the  dispersion  curves  are  virtually  parallel,  with  the  slope  u/k  =  c  2 ,  as  they  tend  to  the 
limiting  velocity  v19  of  shear  waves  in  the  inner  material  as  kh  ->  «.  Figures  2,  3,  4, 

5,  each  contain  two  sets  of  dispersion  curves,  set  (a)  corresponding  to  symmetric  modes 
and  set  (b)  to  antisymmetric  modes  at  angles  60°,  45°,  30°  and  0°  respectively.  As  in 
Figure  1 ,  all  the  branches  except  the  fundamental  antisymmetric  mode  (curves  (b), 
exhibit  a  cut-off  frequency  in  the  long  wavelength  limit  kh  -»  0.  The  curves  again 
display  the  ghost  lines  through  not  to  such  a  marked  extent  as  in  Figure  1  since  the 

number  of  curves  is  now  either  9  or  10  in  each  figure.  Despite  this,  it  is  still  possible 

to  identify  three  ghost  lines  in  Figures  2  and  4  as  opposed  to  two  in  Figures  1 ,  3  and  5 
for  each  of  which  two  of  the  4  speeds  vld>  v  Jd>  vts>  v2S,  become  equal. 

Figures  6,  7  and  8  show  graphs  of  the  group  velocity  Cg,  against  kh.  These  have 

been  obtained  by  numerical  differentiation  of  the  dispersion  curves,  using  a  central 

difference  formula.  Figure  6  corresponds  to  the  First  20  of  the  26  dispersion  curves,  at 

7  =  90°,  shown  in  Figure  1.  On  the  scale  of  the  figure,  the  two  ghost  velocities  are 
vld  =  c,  =  3.25  and  v2d  =  v2J  =  c3  =  2.38  with  the  limiting  velocity  at 
vis  =  c2  =  1.57.  The  curves  clearly  show  the  plateau  and  step  regions  associated  with 
both  ghost  velocities,  displaying  local  maxima  which  approach  c , ,  followed  by  local 
maxima  close  to  c3  and  tending  lo  c,  as  kh  increases.  Figure  7a  shows  the  first  five 
harmonics  of  the  symmetric  modes  at  7  =  60°  and  Figure  7b  gives  the  First  five 

harmonics  of  the  antisymmetric  mode  for  the  same  angle.  At  this  angle  of  propagation 

we  have  v,d  =  3.06,  vjd  =  2.62,  v2J  =  2.20  and  v,,  =  1.81.  None  of  our  curves 
approaches  the  value  of  vld  which  is  not  surprising  since  it  may  be  seen  from  Figure  2 
that  the  first  ghost  line  is  associated  with  harmonics  of  order  greater  than  five.  It  may 
be  seen  that  there  is  one  of  the  five  curves  of  Figure  7a  which  has  a  maximum  at 
Cg  =  2.4,  this  being  the  highest  of  the  five  harmonics  which  may  be  seen  from  2a  to 
have  a  plateau  close  to  the  second  ghost  line.  Both  Figures  7a  and  7b  possess 
harmonics  with  maxima  close  to  Cg  =  2.20,  corresponding  to  the  third  ghost  line 
associated  with  this  angle  of  propagation.  The  curves  shown  in  Figure  8  relate  to  the 
angle  7=0,  the  set  (a)  being  the  first  5  symmetric  modes  and  the  set  (b)  the  first  5 
antisymmetric  modes.  Here,  as  for  7  =  90°,  there  are  only  two  ghost  lines  at 
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v2d  =  c,  =  3.25  and  v,d  =  vls  =  c3  =  2.38.  with  the  lowest  shear  velocity  being 
v!s  .  c,  =  1.57,  in  the  outer  material.  None  of  the  curves  in  8(a)  and  8(b)  have 
maxima  approaching  the  first  ghost  velocity  but  the  maxima  at  the  second  are  very 
clearly  exhibited  in  both  sets.  We  have  remarked  in  Section  3  that  the  local  maxima 
and  minima  of  the  group  velocity  curves  correspond  to  the  wavefronts,  which  are 
expected  to  dominate  the  disturbance  at  large  times.  The  higher  harmonics  shown  in 
Figures  6-8  exhibit  long  flat  maxima  at  the  ghost  velocities,  which  do  not  appear  in  the 
fundamental  modes.  It  is  pointed  out  by  Jones  [9]  for  an  isotropic  plate,  that  the 
residues  associated  with  the  higher  harmonics  are  small  compared  with  those  arising  from 
the  lowest  modes  and  our  numerical  results  show  this  to  be  the  case  here  also. 
Nevertheless,  the  total  contribution  arising  from  the  succession  of  plateau  regions 
associated  with  these  high  harmonics  at  the  ghost  velocities,  will  give  rise  to  precursor 
waves  which  are  not  exhibited  by  the  solution  from  the  fundamental  mode  alone  nor 
from  approximate  plate  theories  designed  to  reproduce  the  fundamental  mode. 

The  remaining  figures  (9-16)  display  the  variation  of  displacements  and  stresses  as 
functions  of  distance  from  the  impact  point  at  different  values  of  the  reduced  time 
T  =  c,t//h  for  propagation  at  an  angle  y  =  60°.  Figures  9a  and  9b  show  the  normal 
displacement  u,  at  the  upper  and  lower  surfaces  respectively.  Each  shows  the  variation 
with  distance  at  times  T  =  10,  20,  30,  40,  50.  Figures  10a  and  10b  show  the  upper 
and  lower  normal  surface  displacements  at  considerably  larger  times,  namely  T  =  100. 
200  and  500.  Figures  11-14,  inclusive,  display  the  normal  displacement  u,,  tangential 
displacement  u2,  the  discontinuous  shear  traction  t,  3  along  the  fibre  direction  and  the 
in-plane  shear  stress  component  t23  respectively,  all  at  time  T  =  40.  Each  figure 
presents  results  at  both  the  upper  and  lower  surfaces,  derived  from  using  the  transforms 
given  in  equation  (43)  at  the  upper  surface  and  their  equivalents  at  the  lower  surface. 

It  may  be  seen  from  equations  (43)  that  the  in-plane  stress  components  1 2  2  and  1 3  3 

have  terms  in  the  transforms  which  involve  kV^  and  these  terms  may  be  evaluated 
directly  from  the  stress  term  t23  by  an  appropriate  scaling.  The  normal  displacements, 
shown  in  Figure  11  are  obtained  by  inverting  the  transform  given  by  equation  (33)  using 
the  Hamming  window.  The  results  show  little  trace  of  a  spurious  oscillation  and  the 
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windowing  appears  to  be  successful.  The  discontinuous  tangential  stress  t,  3  shown  in 
Figure  13  is  in  effect  derived  from  this  displacement  by  differentiation  with  respect 
to  x.  This  serves  to  roughen  the  numerical  results  and  there  is  now  some  indication  of 
a  superimposed  oscillation.  The  tangential  displacement  u2  is  obtained  by  inverting  the 
first  expression  in  equation  (43)  which  involves  multiplying  the  normal  displacement 
transform  by  the  factor  p43/p44.  The  results  shown  in  Figure  12  exhibit  a  considerable 
windowing  effect,  despite  the  use  of  the  Hamming  technique  and  this  is  further 
accentuated  in  Figure  14  since  the  in-plane  tangential  stress  is  essentially  the  x  derivative 
of  the  tangential  displacement. 

In  order  to  examine  the  effects  of  the  harmonics,  we  have  calculated  the  normal 
and  tangential  components  of  displacements  associated  with  the  fundamental  mode  of  the 
symmetric  and  antisymmetric  waves  at  T  =40.  These  are  shown  in  Figure  15  and  it 
may  be  seen  that  each  gives  a  disturbance  commencing  at  x  =  27,  corresponding  to  a 
speed  v  =  2.19  in  the  units  we  are  employing.  These  curves  may  be  compared  with 
those  sho^h  in  Figure  11  wich  are  calculated  using  a  total  of  18  modes  (9  each  of  the 
symmetric  and  antisymmetric)  and  which  show  a  disturbance  commencing  at  x  =  32 
corresponding  to  v  =  2.60.  Finally,  in  Figure  16  we  show  the  normal  and  tangential 
upper  surface  displacements  at  T  =40°  calculated  using  10  modes  only  (5  each  of  the 
symmetric  and  antisymmetric).  We  note  that  these  produce  a  disturbance  at  only  slightly 
larger  values  of  x  than  the  two  fundamental  modes.  This  is  consistent  with  the  fact  that 
the  highest  ghost  speed  is  only  displayed  by  the  harmonics  of  order  higher  than  five  in 
each  mode  of  wave  motion.  We  further  note  that  these  curves  show  little  evidence  of 
oscillation  in  the  tangential  displacement  as  compared  with  that  displayed  in  Figure  12. 
This  is  consistent  with  the  fact  that  all  ten  first  harmonics  have  crossed  the  lowest  ghost 
line  before  the  cut  off  wave  number  k  is  reached.  The  motion  associated  with  each  of 
these  harmonics  is  then  a  propagating  shear  disturbance  largely  confined  to  the  inner 
core  and  these  harmonics  make  no  contribution  to  the  surface  displacements  at  values 
close  to  k  so  that  there  is  no  effect  to  be  expected  from  curtailing  the  integrals. 
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Dispersion  curves  for  propagation  at  y  =  90°. 


Figure  2b.  Dispersion  curves  for  antisymmetric  modes  of  propagation  at  7  =  60° 
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Dispersion  curves  for  symmetric  modes  of  propagation  at  y  =  0° 
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Figure  8b  Group  velocity  curves  for  first  5  antisymmetric  modes  of  propagation 

at  7  =  0° 
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Figure  15  (a)  Normal  displacement  (b)  Tangential  displacement 


at  the  upper  surface  arising  from  fundamental  modes  of 
(i)  symmetric  and  (ii)  antisymmetric  motion  at  T  =  40. 


Figure  18  Geometry  of  impact  problem 
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